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Abstract. Iterative algorithms aimed at solving some problems are discussed. 

UFor certain problems, such as finding a common point in the intersection of a finite 
number of convex sets, there often exist iterative algorithms that impose very little 
Q ~) , demand on computer resources. For other problems, such as finding that point in 

• i the intersection at which the value of a given function is optimal, algorithms tend 

^_ to need more computer memory and longer execution time. A methodology is 

presented whose aim is to produce automatically for an iterative algorithm of the 
first kind a "superiorized version" of it that retains its computational efficiency 
but nevertheless goes a long way towards solving an optimization problem. This 
is possible to do if the original algorithm is "perturbation resilient," which is 
shown to be the case for various projection algorithms for solving the consistent 
convex feasibility problem. The superiorized versions of such algorithms use 
^ perturbations that drive the process in the direction of the optimizer of the given 

^^\ function. After presenting these intuitive ideas in a precise mathematical form, 

VP , they are illustrated in image reconstruction from projections for two different 

f^*^ ■ projection algorithms superiorized for the function whose value is the total 

f^ ■ variation of the image. 

in 

Keywords: iterative algorithms, convex feasibility problem, superiorization, 
perturbation resilience, projection methods 



o 



- 1—1 

X 



AMS classification scheme numbers: 65Y20, 68W25, 90C06, 90C25, 68U10 



Submitted to: Inverse Problems 



Perturbation Resilience and Superiorization 2 

1. Introduction 

We first motivate and describe our ideas in a not fully general context, in which 
superiorization is envisioned as lying in-between the methodologies of optimization 
and of feasibility seeking. With a feasible solution one settles for a point that just 
fulfills a set of constraints, whereas solving a constrained optimization problem calls for 
finding a feasible point that optimizes a given objective function. Generally speaking, 
optimization is logically and computationally a more demanding task than that of 
finding just any feasible point. We show that, without employing an optimization 
algorithm, it is possible to use certain iterative methods, designed for (the less 
demanding) feasibility problems, in a way that will steer the iterates toward a point 
that is superior, but not necessarily optimal, in a well-defined sense. The advantage 
of superiorization is that it allows us to solve significant problems by using powerful 
feasibility seeking methods, see, e.g., |10] and references therein, and reach a superior 
feasible point without resorting to optimization techniques. We now explain this with 
more details. 

Many significant real-world problems are modeled by constraints that force the 
sought-after solution point to fulfill conditions imposed by the physical nature of the 
problem. Such a modeling approach often leads to a convex feasibility problem of the 
form 

I 
find x* e C = p| C it (1) 

8=1 

where the sets C, C R J are closed convex subsets of the Euclidean space K J , see 
[H OH OEH] or P21 Chapter 5] for this broad topic. In many real-world problems the 
underlying system is very large (huge values of I and J) and often very sparse. In 
these circumstances projection methods have proved to be effective. They are iterative 
algorithms that use projections onto sets while relying on the general principle that 
when a family of closed and convex sets is present, then projections onto the individual 
sets are easier to perform than projections onto other sets, such as their intersection 
as in ([1]), that are derived from them. 

Projection methods can have various algorithmic structures (some of which 
are particularly suitable for parallel computing) and they also possess desirable 
convergence properties and good initial behavior patterns [21 H3 El [TSl HH [2Z1 132] • 
The main advantage of projection methods, which makes them successful in real- 
world applications, is computational. They commonly have the ability to handle 
huge-size problems of dimensions beyond which more sophisticated methods cease 
to be efficient or even applicable due to memory requirements. (For a justification 
of this claim see the various examples provided in [TU].) This is so because the 
building bricks of a projection algorithm (which are the projections onto the given 
individual sets) are easy to perform, and because the algorithmic structure is either 
sequential or simultaneous, or in-between, as in the block-iterative projection methods 
or in the more recently invented string-averaging projection methods. The number 
of sets used simultaneously in each iteration in block-iterative methods and the 
number and lengths of strings used in each iteration in string-averaging methods 
are variable, which provides great flexibility in matching the implementation of the 
algorithm with the parallel architecture at hand; for block- iterative methods see, e.g., 
[U 02 [H OH OH 022 OH 0E2 US 02] and for string-averaging methods see, e.g., 

[iigmnaiHiiaosii!!]. 
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The key to superiorization is our recent discovery [HI 1231 128| that two principal 
prototypical algorithmic schemes of projection methods: string-averaging projections 
(SAP) and block- iterative projections (BIP), which include as special cases a variety 
of projection methods for the convex feasibility problem, are bounded perturbations 
resilient in the sense that the convergence of sequences generated by them continues 
to hold even if the iterates are perturbed in every iteration. We harness this resilience 
to bounded perturbations to steer the iterates to not just any feasible point but to a 
superior (in a well-defined sense) feasible point of ([lj. 

Our motivation is the desire to create a new methodology that will significantly 
improve methods for the solution of inverse problems in image reconstruction from 
projections, intensity-modulated radiation/proton therapy (IMRT/IMPT) and in 
other real- world problems such as electron microscopy (EM). Our work [HI H3], as 
well as the examples given below, indicate that our objective is achievable and show 
how algorithms can incorporate perturbations in order to perform superiorization. 

The superiorization methodology has in fact broader applicability than what has 
been discussed until now and its mathematical specification in the next section reflects 
this. However, all our specific examples will be chosen from the field that we used as 
our motivation in this introductory section. 

2. Specification of the superiorization methodology 

The superiorization principle relies on the bounded perturbation resilience of 
algorithms. Therefore we define this notion next in a general setting within R J . 

We introduce the notion of a problem structure (T, Vr), where T is a nonempty 
problem set and Vr is a function on T such that, for all T £ T, Vrr ■ R J — > R+, 
where R + is the set of nonnegative real numbers. Intuitively we think of Vtt (x) as a 
measure of how "far" x is from being a solution of T. In fact, we call x a solution of 
T if Vr T (x) = 0. 

For example, for the convex feasibility problem ([1]) 

T = {{Ci, . . . , Ci}\ I is a positive integer and, for 1 < i < I, , . 
Ci is a closed convex subset of IR 1 - 7 } 

and 



\ 



•£(d(x,C,)) 2 , (3) 



where d(x,Ci) is the Euclidean distance of x from the set Cj. Clearly, in this case 
a; is a solution of {Ci, . . . , C/} as defined in the previous paragraph if, and only if, 
xeCas defined in (TIJ. 

Definition 1. An algorithm P for (T,Vr) assigns to each T e Tan algorithmic 
operator Pt : R' 7 — > K' 7 - P is said to be bounded perturbations resilient if, for all 

(Pt) x } converges to a solution 

J fe=0 

of T for all x £ R J , then any sequence {a; fc }, of points in R J also converges to a 
solution of T provided that, for all k > 0, 

x k+1 = P T {x k + p k v k ) , (4) 



Perturbation Resilience and Superiorization 4 

where /3kV k are bounded perturbations, meaning that /3k are real nonnegative numbers 

oo 

such that 2_j@ k < °° an< ^ * ne sequence {i'' c },,_ r) is bounded. 

fc=0 

We give next specific instances of bounded perturbations resilient algorithms for 
solving the convex feasibility problem as in © and ([3]), from the classes of SAP and 
BIP methods. We do this by defining P{c 1 ....,c T } f° r an arbitrary but fixed element 
{Ci,...,C/} of T of © for the different algorithms P. For any nonempty closed 
convex subset M of R J and any x 6 M. J , the orthogonal projection of x onto M is the 
point in M that is nearest (by the Euclidean distance) to x; it is denoted by PmX. 

To define P{Ci ...,Ci} f° r the SAP instances, we make use of index vectors, which 
are nonempty ordered sets t = {t\, . . . , ijy), where N is an arbitrary positive integer, 
whose elements t n are in the set {1, ...,/} . For an index vector t we define the composite 
operator 

P[t]=P CtN ---Pc H . (5) 

A finite set fi of index vectors is called fit if, for each i G {1,...,/}, there exists 
t = (t\, . . . , ijv) £ £1 such that t n = i for some n 6 {1, ..., N} . If tt is a fit set of 
index vectors, then a function lu : O — > IR++ = (0, oo) is called a fit weight function if 
Stefi w (*) = 1- A pair (f2, uS) consisting of a fit set of index vectors and a fit weight 
function defined on it was called an amalgamator in [B]. For each amalgamator (f2, w) , 
we define the algorithmic operator P{Ci,...,C7} : K J — >• M J by 

p { c 1 ,..,c f } a; =E w ( t ) p W a; - ( 6 ) 

ten 
For this algorithmic operator we have the following bounded perturbations resilience 
theorem. 

Theorem 1. [6, Section II] If C of flj) is nonempty, {/3fc}^L is a sequence of non- 
negative real numbers such that X)fc°=o^ fc *"• °° an( ^ \ v ih—n * s a bounded sequence 
of points in M. J , then for any amalgamator (£l,ui) and any x° £ M. J , the sequence 
{x }i_ generated by 



x 



fe+i _ 



P {Cl ,...^}( a;fe + / 3 ^ fe ) I Vfc >°> ( 7 ) 

converges, and its limit is in C. (The statement of this theorem in [5J is for positive 
/?fcS, but the proof given there applies to nonnegative /3feS.) 

Corollary 1. For any amalgamator ($7,w), the algorithm P defined by the algo- 
rithmic operator Pr ClCl y is bounded perturbations resilient. 

f k "I °° 

Proof. Assume that for T = {Ci, . . . , C/} the sequence < (Pt) a? f converges to 

a solution of T for all x G M J . This implies, in particular, that C of (pj is nonempty. 
By Definition 1, we need to show that any sequence {a; fc }, of points in K. J also 
converges to a solution of T provided that, for all k > 0, ^ is satisfied when the (3kV k 
are bounded perturbations. Under our assumptions, this follows from Theorem 1. □ 

Next we look at a member of the family of BIP methods. Considering the convex 
feasibility problem ([lj, for 1 < u < U, let B u be an ordered set (b u ,i, ■ ■ ■ ,&«,|b u |) of 
elements of {1, . . . , 1} (\B U \ denotes the cardinality of B u ). We call such a B u a block 
and define the (composite) algorithmic operator Q{Ci,...,Cj} : K J — > K J by 

Q{Ci,...,Cj} = Qu ■ ••Qi) (8) 
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where, for x £ K J and 1 < u < U , 



Q^^E^ + ^IT^' ( 9 ) 



and 



i? = max{|£ u | | 1 <u< U} . (10) 

The iterative procedure x k+1 — Q{d,...,Ci} x 1S a member of the family of BIP 
methods. For this algorithmic operator we have the following bounded perturbations 
resilience theorem. 

Theorem 2. |23| If C of ^ is nonempty, {1,...,/} = \J U=1 B U , {/3k}^ =0 is a 
sequence of nonnegative real numbers such that X)fc°=o @k < °°i \ v it— n ^ e a bounded 
sequence of points in M. J , then for any x° £ M. J , the sequence {x k \ generated by 

x k+1 = Q { c 1 ,....c,} {x k + hv k ) , Vfc > 0, (11) 

converges, and its limit is in C. (This is a special case of Theorem 2 in |23j given here 
without a relaxation parameter. Also, that theorem is stated for positive /3fcS, but the 
proof given there applies to nonnegative /3feS.) 

Corollary 2. The algorithm Q defined by the algorithmic operator Q{Ci,...,Cr} * s 

bounded perturbations resilient. 

Proof. Replace in the proof of Corollary 1 P by Q and Theorem 1 by Theorem 2. □ 

Further bounded perturbations resilience theorems are available in a Banach 
space setting, see pj [8]. Thus the theory of bounded perturbations resilient 
algorithms already contains some solid mathematical results. As opposed to this, 
the superiorization theory that we present next is at the stage of being a collection 
of heuristic ideas, a full mathematical theory still needs to be developed. However, 
there are practical demonstrations of its potential usefulness; see [51 [231 HH] and the 
illustrations in Section [3] below. 

For a problem structure (T,Vr), T £ T, e £ M ++ and a sequence S = {x k } k=Q 
of points in M J , we use 0(T 7 e 7 S) to denote the x £ R J that has the the following 
properties: Vrxix) < e and there is a nonnegative integer K such that x K = x 
and, for all nonnegative integers £ < K, Vrx (as ) > e. Clearly, if there is such an 
x, then it is unique. If there is no such x, then we say that O (T, e, S) is undefined. 
The intuition behind this definition is the following: if we think of S as the (infinite) 
sequence of points that is produced by an algorithm (intended for the problem T) 
without a termination criterion, then O (T, e, S) is the output produced by that al- 
gorithm when we add to it instructions that makes it terminate as soon as it reaches 
a point at which the value oiVrx is not greater than e. The following result is obvious. 

Lemma 1. If Vtt is continuous and the sequence S converges to a solution of T , 
then O (T, e, S) is defined and Vrx (O (T, e, S)) < e. 

Given an algorithm P for a problem structure (T, Vr), a T £ T and an 

x £ R J , let R(T,x) = \(P T ) k x\ . For a function <f> : R J -4 R, the supe- 

L J k=o 

riorization methodology should provide us with an algorithm that produces a se- 
quence S(T,x,(j)) — {a; fe } fc such that for any e £ M ++ and x £ R J for which 
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Pr T (x) > e and 0(T,s,R(T,x)) is defined, O (T,e, S (T, x, </>)) is also defined and 
cj> (O (T, e, S (T, S, (f>))) < cj)(0 (T, e, E (T, x))). This is of course too ambitious in its 
full generality and so here we analyze only a special case, but one that is still quite 
general. We now list our assumptions for the special case for which we discuss details 
of the superiorization methodology. 

Assumptions 

(i) (T, Vr) is a problem structure such that Vtt is continuous for all T G T. 
(ii) P is a bounded perturbation resilient algorithm for (T, Vr) such that, for all T G 
T. Pt is continuous and, if x is not a solution of T, then Vrx (Ptx)) < Vrr (x). 
(iii) <j) is a convex function. 

We now describe, under these assumptions, the algorithm to produce the sequence 
S{T,x,<j>) = {x k }™ =Q . 

The algorithm assumes that we have available a summable sequence {"fe}g = Q 
of positive real numbers. It is easy to generate such sequences; e.g., we can use 
jg = a e , where < a < 1. The algorithm generates, simultaneously with the sequence 
{x },_ n , sequences {«},_„ and {/Sfcl^Lg- ^ ne l & tter will be generated as a subse- 
quence of {-fe }^ . Clearly, the resulting sequence {/3fe}^L of positive real numbers 
will be summable. We first specify the algorithm and then discuss it. The algorithm 
depends on the specified x, </>, {7f}^ 0! Vrx and Py. It makes use of a logical variable 
called continue and also of the concept of a subgradient of the convex function <f>. 

Superiorized Version of Algorithm P 



(i 
(ii 
(iii 

(iv 
(v 

(vi 

(vii 
(viii 

(ix 
(x 

(xi 

(xii 
(xiii 
(xiv 

(xv 
(xvi 
(xvii 



set k = 
set x = x 
set £ = 
repeat 

set g to a subgradient of <j> at x k 

if llffll>0 

then set v = —g/ \\g\\ 

else set v k = g 
set continue = true 
while continue 

set @ k = 7£ 

set y = x k + PkV k 

if 4>{y) <(f> (x k ) and Vr T (Pry) < P?"T (^' £ ) then 
set a; fc+1 = P T y 
set continue = false 

set ^ = £ + 1 
set fc = k + 1 
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Sometimes it is useful to emphasize the function <f> for which we are superiorizing, in 
which case we refer to the algorithm above as the <f>- super iorized version of algorithm 
P. It is important to bear in mind that the sequence S produced by the algorithm 
depends also on the initial point x, the selection of the subgradient in Line (v) of 
the algorithm, the summable sequence {7^}^ =0 , and the problem T. In addition, the 
output O (T, s, S) of the algorithm depends on the stopping criterion e. 

Theorem 3. Under the Assumptions listed above, the Superiorized Version of Al- 
gorithm P will produce a sequence S (T, x,(f>) of points in M. J that either contains a 

solution ofT or is infinite. In the latter case, if the sequence \ (Pt) x \ converges 

L J fc=o 

to a solution ofT for all x £ R" 7 , then, for any e £ R++, O (T, e, S (T, x, </>)) is defined 

and 4>{0 (T,e,S (T,x,4>))) < e. 

Proof. Assume that the sequence S (T, x, 4>) produced by the Superiorized Version 
of Algorithm P dos not contain a solution of T. We first show that in this case the 
algorithm generates an infinite sequence {x k } . This is equivalent to saying that, 
for any x k that has been generated already, the condition in Line (xiii) of the al- 
gorithm will be satisfied sooner or later (and hence x k+1 will be generated). This 
needs to happen, because as long as the condition is not satisfied we keep reset- 
ting (in Line (xi)) the value of /3fe to 7^, with ever increasing values of £. However, 
{li}i=o ^ s a summable sequence of positive real numbers, and so 7^ is guaranteed 
to be arbitrarily small if £ is sufficiently large. Since v k is either a unit vector in 
the direction of the negative subgradient of the convex function <j> at x k or is the 
zero vector (see Lines (v)— (viii)), (j) ( x + Pk v ) < <f> ( x ) must be satisfied if the 
positive number /3fc is small enough. Also, since Vtt (Pt^^) < Vtt (x k ) and Pt 
and Vtt are continuous (Assumptions (ii) and (i), respectively), we also have that 
Vrr (Pr (x k + fikV k )) < Vrx (x k ) if j3k is small enough. This completes the proof 
that the condition in Line (xiii) of the algorithm will be satisfied and so the algorithm 
will generate an infinite sequence S (T, x, cj>). Observing that we have already demon- 
strated that the j3kV k are bounded perturbations, and comparing (|4]) with Lines (xii) 

and (xiv), we see that (by the bounded perturbation resilience of P) the assumption 

{, > 00 
(Pt) x > converges to a solution of T for all x G M J implies 
i k=0 

that S (T, x, <p)) also converges to a solution of T . Thus, applying Lemma 1 we obtain 
the final claim of the theorem. □ 

Unfortunately, this theorem does not go far enough. To demonstrate that a 
methodology leads to superiorization we should be proving (under some assumptions) 
a result like <f> (O (T, e, S (T, x, (/>))) < <j)(0 (T, e, R (T, x))) in place of the weaker result 
at the end of the statement of the theorem. Currently we do not have any such proofs 
and so we are restricted to providing practical demonstrations that our methodology 
leads to superiorization in the desired sense. In the next section we provide such 
demonstrations for the Superiorized Version of Algorithm P, for two different Ps. 

3. Illustrations of the superiorization methodology 

We illustrate the superiorization methodology on a problem of reconstructing a head 
cross-section (based on Figure 4.6(a) of [27]) from its projections using both an SAP 
and a BIP algorithm. (All the computational work reported in this section was done 
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Table 1: Values of TV for the outputs of the various algorithms. The second column 
is for the superiorized versions and the third column is for the original versions. 

Algorithm <f>(0(T,e,S(T,x,(f>))) (O (T, e, R (T, x)))) 

Variant of ART 441.50 1,296.44 

Variant of BIP 444.15 1,286.44 



using SNARK09 [24]; the phantom, the data, the reconstructions and displays were 
all generated within this same framework.) Figure QJa) shows a 243 x 243 digitization 
of the head phantom with J = 59, 049 pixels. An x £ M. J is interpreted as a vector 
of pixel values, whose components represent the average X-ray linear attenuation 
coefficients (measured per centimeter) within the 59, 049 pixels. Each pixel is of size 
0.0752 x 0.0752 (measured in centimeters). The pixel values range from to 0.5639. 
For display purposes, any value below 0.204 is shown as black (gray value 0) and any 
value above 0.21675 is shown as white (gray value 255), with a linear mapping of the 
pixel values into gray values in between (the same convention is used in displaying 
reconstructed images in Figures QJb)-(e)). 

Data were collected by calculating line integrals across the digitized image for 82 
sets of equally spaced parallel lines, with / = 25,452 lines in total. Each data item 
determines a hyperplane in WL J . Since the digitized phantom lies in the intersection of 
all the hyperplanes, we have here an instance of the convex feasibility problem with a 
nonempty C, satisfying the first condition of the statements of Theorems 1 and 2. 

For our illustration, we chose the SAP algorithm P{Ci,....Cr} as determined by 
(JSJl-jl]) with fl — {(1, ... , I)} and u (1, . . . , I) = 1. This is a classical method that in 
tomography would be considered a variant of the algebraic reconstruction techniques 
(ART) [S7J Chapter 11]. For the BIP algorithm we chose Q{Ci,...,Ci} as determined 
by ([5])- (fTU]) with U — 82 and each block corresponding to one of the 82 sets of parallel 
lines along which the data are collected. 

The function <fi for which we superiorized is defined so that, for any x £ M. J , (f> (x) 
is the total variation (TV) of the corresponding 243 x 243 image. If the pixel values 
of this image are q g ,h, then the value of the TV is defined to be 



Yl 5Z y(l9+hh - q g ,hf + {q g ,h+i - qgjif- (12) 

g=lh=l 

For the TV-Superiorized Versions of the Algorithms P{Ci,...,cM an d Q{Ci,...,Ci} of 
the previous paragraph we selected x to be the origin (the vector of all zeros) and 
7£ = 0.999 f . Also , we set e — 0.01 for the stopping criterion, which is small compared 
to the Vr T of the initial point (Vr T (x) = 330.208). 

For each of the four algorithms (P{Ci,...,Ci\: Q{Ci,...,c7r} an d their TV- 
superiorized versions), the sequence S that is produced by it is such that the output 
0(T,e,S) is defined; see Figures HJb)-(e) for the images that correspond to these 
outputs. Clearly, the superiorized reconstructions in Figures QJc) an d ( e ) ar e visually 
superior to their not superiorized versions in Figure [Tfb) and (d) , respectively. More 
importantly from the point of view of our theory, consider TableQ] As stated in the last 
paragraph of the previous section, we would like to have that </> (O (T, e, S (T, x,(f>))) < 
4> (O (T, e, R (T, x))). While we are not able to prove that this is the case in general, 
Table [1] clearly shows it to be the case for the two algorithms discussed in this section. 
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Figure 1: A head phantom (a) and its reconstructions from underdetermined 
consistent data obtained for 82 views using: (b) a variant of ART, (c) TV-superiorized 
version of the same variant of ART, (d) a block-iterative projection method, and (e) 
TV-superiorized version of the same block-iterative projection method. The same 
initial point and stopping criterion were used in all cases; see the text for details. 
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A final important point that is illustrated by the experiments in this section 
is that, from the practical point of view, TV-superiorization is as useful as TV- 
optimization. This is because a realistic phantom, such as the one in Figure 
[lja), is unlikely to be TV- minimizing subject to the constraints provided by the 
measurements. In fact, the TV value of our phantom is 450.53, which is larger than 
that for either of the TV-superiorized reconstructions in the second column of Table 
[T] While an optimization method should be able to find an image with a lower TV 
value, there is no practical point for doing that. The underlying aim of what we 
are doing is to estimate the phantom from the data and producing an image whose 
TV value is further from the TV value of the phantom than that of our superiorized 
reconstructions is unlikely to be helpful towards achieving this aim. 

4. Discussion and conclusions 

Stability of algorithms under perturbations is generally studied in numerical analysis 
with the aim of proving that an algorithm is stable so that it can "endure" all kinds of 
imperfections in the data or in the computational performance. Here we have taken 
a proactive approach designed to extract specific benefits from the kind of stability 
that we term perturbation resilience. We have been able to do this in a context that 
includes, but is much more general than, feasibility-optimization for intersections of 
convex sets. 

Our premise has been that (1) there is available a bounded perturbations resilient 
iterative algorithm that solves efficiently certain type of problems and (2) we desire 
to make use of perturbations to find for these problems solutions that, according to 
some criterion, are superior to the ones to which we would get without employing 
perturbations. To accomplish this one must have a way of introducing perturbations 
that take into account the criterion according to which we wish to "superiorize" the 
solutions of the problems. 

We have set forth the fundamental principle, have given some mathematical 
formulations and results, and have shown potential benefits (in the field of image 
reconstruction from projections). However, the superiorization methodology needs to 
be studied further from the mathematical, algorithmic and computational points of 
view in order to unveil its general applicability to inverse problems. As algorithms 
are developed and tested a dialog on algorithmic developments must be accompanied 
by mathematical validation and applications to simulated and real data from various 
relevant fields of applications. 

Validating the concept means proving precise statements about the behavior of 
iterates {^ fe },._ n generated by the superiorized versions of algorithms. Under what 
conditions do they converge? Can their limit points be characterized? How would 
different choices of the perturbation coefficients /?& and the perturbation vectors v k 
affect the superiorization process? Can different schemes for generating the /3fcS be 
developed, implemented, investigated? Enlarging the arsenal of bounded perturbation 
resilience algorithms means generalizing existing proofs of bounded perturbations 
resiliency of algorithms and developing new theories that will bring more algorithms 
into the family of bounded perturbations resilient algorithms. Further developments 
should include the problem of finding a common fixed point of a family of operators 
(a direct generalization of the convex feasibility problem) and studying the behavior 
of superiorization algorithms in inconsistent situations when the underlying solution 
set is empty. Thus we view the material in this paper as only an initial step in a 
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promising new field of endeavor for solving inverse problems. 
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